require(Rcpp)
require(RcppArmadillo)
require(coda)
require(lmtest)
require(sandwich)
require(RColorBrewer)
Set1 <- brewer.pal(9, "Set1")
windowsFonts(Helvetica = windowsFont("Helvetica"))

attribute.label <- c("Parents", "Identity", "Activity", 
                     "Policy", "Habit", "Intention")
level.label <- list(c("Independent", "Same party", "Opponent party"), 
                    c("Does not matter", "Personal insult"), 
                    c("Never", "Rarely", "Sometimes", "Often", "Very often"), 
                    c("Does not approve", "Partially approves", "Completely approves"), 
                    c("Always abstained", "Sometimes abstained", 
                      "Sometimes deviated", "Always loyal"), 
                    c("Abstain", "Same party", "Opponent party"))

# Gibbs sampler code for a finite mixture linear regression model
# with the multinomial logit (MNL) classification
sourceCpp("mixture_reg_with_MNL.cpp")

# Gibbs sampler code for a finite mixture linear regression model
# without the MNL classification
sourceCpp("mixture_reg_wo_MNL.cpp")

# function to create a mcmc.list object
mcmc.list.create <- function(mcmc.result, parameter, label = NULL, 
                             group = NULL, burnin = 0) {
  combined.list <- mcmc.list()
  if (burnin > 0) {
    sampled <- c((burnin + 1):dim(mcmc.result[[1]][[1]])[3])
  } else {
    sampled <- c(1:dim(mcmc.result[[1]][[1]])[3])
  }
  if (parameter %in% c("beta", "beta_r", "beta_c")) {
    if (is.null(label)) {
      for (i in 1:length(mcmc.result)) {
        combined.list[[i]] <- mcmc(t(mcmc.result[[i]][[parameter]][, , sampled]))
      }
    } else {
      for (i in 1:length(mcmc.result)) {
        combined.list[[i]] <- mcmc(t(mcmc.result[[i]][[parameter]][, label[i, group], sampled]))
      }
    }
  } else {
    if (is.null(label)) {
      for (i in 1:length(mcmc.result)) {
        combined.list[[i]] <- mcmc(mcmc.result[[i]][[parameter]][sampled, ])
      }
    } else {
      for (i in 1:length(mcmc.result)) {
        combined.list[[i]] <- mcmc(mcmc.result[[i]][[parameter]][sampled, label[i, group]])
      }
    }
  }
  combined.list 
}

# function to relabel the MCMC chains of parameters other than gamma
labeling <- function(mcmc.result, burnin = 0) {
  MCMCsample <- list()
  for (i in 1:length(mcmc.result)) {
    if (burnin > 0) {
      MCMCsample[[i]] <- mcmc.result[[i]]$G[, , -c(1:burnin)]
    } else {
      MCMCsample[[i]] <- mcmc.result[[i]]$G
    }
  }
  l <- length(MCMCsample)
  g <- dim(MCMCsample[[1]])[2]
  label <- matrix(NA, l, g)
  for (i in 1:l) {
    group.prob <- apply(MCMCsample[[i]], 2, mean)
    for (j in 1:g) {
      label[i, j] <- which(c(g - rank(group.prob) + 1) == j)
    }
  }
  label
}

# function to relabel the MCMC chains of gamma
gamma.labeling <- function(mcmc.result, label, burnin = 0) {
  MCMCsample <- list()
  for (i in 1:length(mcmc.result)) {
    if (burnin > 0) {
      MCMCsample[[i]] <- mcmc.result[[i]]$gamma[, , -c(1:burnin)]
    } else {
      MCMCsample[[i]] <- mcmc.result[[i]]$gamma
    }
  }
  output <- list()
  for (i in 1:nrow(label)) {
    temp.list <- array(NA, dim(MCMCsample[[1]]))
    for (j in 1:dim(MCMCsample[[1]])[1]) {
      if (dim(MCMCsample[[1]])[2] == 1) {
        temp.gamma <- cbind(MCMCsample[[i]][j, , ], 0)
      } else {
        temp.gamma <- cbind(t(MCMCsample[[i]][j, , ]), 0)
      }
      for (k in 1:dim(MCMCsample[[1]])[2]) {
        temp.list[j, k, ] <- exp(temp.gamma[, label[i, k + 1]]) / 
          exp(temp.gamma[, label[i, 1]])
      }
    }
    output[[i]] <- temp.list
  }
  output
}

# function to compute the summary statistics of posterior distributions
posterior.summary <- function(x, CI = TRUE) {
  posterior.mean <- mean(x)
  if (CI == TRUE) {
    HPD.interval <- HPDinterval(as.mcmc(x))
    c(posterior.mean, HPD.interval)
  } else {
    posterior.mean
  }
}

# function to compute the WAIC
WAIC <- function(loglik) {
  lppd <- sum(log(colMeans(exp(loglik))))
  p.WAIC2 <- sum(apply(loglik, 2, var))
  -2 * (lppd - p.WAIC2)
}

# function to conduct postestimation simulation for the MNL
MNL.first.dif <- function(MCMCsample, label, data, colnum, low, high, 
                          dummy, dummycol = NULL, prob) {
  sim.data <- cbind(1, data)
  if (dummy == TRUE) {
    for (i in 1:length(dummycol)) {
      sim.data[, dummycol[i] + 1] <- min(data[, dummycol[i]])
    }
  }
  result <- list()
  for (i in 1:length(MCMCsample)) {
    result[[i]] <- matrix(NA, dim(MCMCsample[[i]])[3], 
                          dim(MCMCsample[[i]])[2] + 1)
    for (j in 1:dim(MCMCsample[[i]])[3]) {
      sim.data[, colnum + 1] <- low
      expzgamma <- exp(sim.data %*% MCMCsample[[i]][, , j])
      prob.low <- cbind(expzgamma, 1) / (1 + rowSums(expzgamma))
      sim.data[, colnum + 1] <- high
      expzgamma <- exp(sim.data %*% MCMCsample[[i]][, , j])
      prob.high <- cbind(expzgamma, 1) / (1 + rowSums(expzgamma))
      result[[i]][j, ] <- colMeans(prob.high - prob.low)
    }
  }
  combined.result <- c()
  for (i in 1:length(MCMCsample)) {
    combined.result <- rbind(combined.result, result[[i]][, label[i, ]])
  }
  cbind(colMeans(combined.result), 
        HPDinterval(as.mcmc(combined.result), prob = prob))
}

# function to compute standard errors
se <- function(x) {
  sd(x) / sqrt(length(x))
}

#### preliminary analysis excluding additional satisficers ####
raw.data <- read.csv("raw_data.csv")
load("attentive_data.Rdata")
load("clean_data.Rdata")

# additional satisficers, who chose the same option in their rating answers
additional.satisficers <- tapply(clean.data$rating, 
                                 clean.data$respondent.id, 
                                 function(x) length(table(x)) == 1)
additional.satisficers.id <- unique(clean.data$respondent.id)[additional.satisficers]

# additional satisficers' responses
for (i in 1:length(additional.satisficers.id)) {
  print(clean.data$rating[clean.data$respondent.id == additional.satisficers.id[i]])
}

# drop those who with missing values and additional satisficers
omit.na.data <- subset(clean.data, 
                       (! is.na(education)) & (! is.na(black)) & 
                         (! respondent.id %in% additional.satisficers.id))

# preliminary estimation of a linear model
preliminary.lm <- lm(rating ~ a.parents + b.identity + c.activity + 
                       d.policy + e.habit + f.intention, 
                     data = omit.na.data, x = TRUE)

# individual mean of rating responses for within-individual standardization
average.rating <- rep(NA, nrow(raw.data))
for (i in unique(omit.na.data$respondent.id)) {
  average.rating[i] <- mean(clean.data$rating[clean.data$respondent.id == i])
}
omit.na.data$average.rating <- average.rating[omit.na.data$respondent.id]

# data objects used to estimate finite mixture linear regression models
y.r <- omit.na.data$rating - omit.na.data$average.rating
y.c <- omit.na.data$choice
x <- preliminary.lm$x[, -1]
# covariates
merged.data <- merge(data.frame(ID = unique(omit.na.data$respondent.id)), 
                     attentive.data, by = "ID")
z <- cbind(merged.data$gender, merged.data$age, merged.data$education, 
           merged.data$black, merged.data$other.race, 
           diag(4)[merged.data$region, -1], diag(3)[merged.data$PID, -3], 
           merged.data$knowledge, merged.data$ideology)
id <- as.numeric(as.factor(omit.na.data$respondent.id))
l <- length(unique(id))

#### estimate finite mixture linear regression models ####
# one-group model
g <- 1
mixture.with.MNL.1.result <- list()
for (i in 1:3) {
  set.seed(100 + 10 * i)
  mixture.with.MNL.1.result[[i]] <- 
    mixture_reg_MNL(y_r = y.r, y_c = y.c, x = x, z = z, 
                    id = id, l = l, g = g, 
                    b0 = matrix(0, ncol(x) + 1, g), 
                    Bvar = matrix(100, ncol(x) + 1, g), 
                    v0 = rep(0.01, g), s0 = rep(0.01, g), 
                    g0 = matrix(rep(0, (ncol(z) + 1) * (g - 1)), ncol(z) + 1, g - 1), 
                    G0 = matrix(rep(100, (ncol(z) + 1) * (g - 1)), ncol(z) + 1, g - 1), 
                    tune = matrix(rep(0.1, (ncol(z) + 1) * (g - 1)), ncol(z) + 1, g - 1), 
                    burnin = 5000, iter = 5000, thin = 50)
}

# save(mixture.with.MNL.1.result, file = "mixture_with_MNL_1_result.Rdata")
# load("mixture_with_MNL_1_result.Rdata")

# convergence check
sum(gelman.diag(mcmc.list.create(mixture.with.MNL.1.result, "beta_r"), 
                multivariate = FALSE)$psrf[, 1] > 1.1)
sum(gelman.diag(mcmc.list.create(mixture.with.MNL.1.result, "beta_c"), 
                multivariate = FALSE)$psrf[, 1] > 1.1)
sum(gelman.diag(mcmc.list.create(mixture.with.MNL.1.result, "sigma_r"), 
                multivariate = FALSE)$psrf[, 1] > 1.1)
sum(gelman.diag(mcmc.list.create(mixture.with.MNL.1.result, "sigma_c"), 
                multivariate = FALSE)$psrf[, 1] > 1.1)

# two-group model
g <- 2
mixture.with.MNL.2.result <- list()
for (i in 1:3) {
  set.seed(200 + 10 * i)
  mixture.with.MNL.2.result[[i]] <- 
    mixture_reg_MNL(y_r = y.r, y_c = y.c, x = x, z = z, 
                    id = id, l = l, g = g, 
                    b0 = matrix(0, ncol(x) + 1, g), 
                    Bvar = matrix(100, ncol(x) + 1, g), 
                    v0 = rep(0.01, g), s0 = rep(0.01, g), 
                    g0 = matrix(rep(0, (ncol(z) + 1) * (g - 1)), ncol(z) + 1, g - 1), 
                    G0 = matrix(rep(100, (ncol(z) + 1) * (g - 1)), ncol(z) + 1, g - 1), 
                    tune = matrix(rep(0.1, (ncol(z) + 1) * (g - 1)), ncol(z) + 1, g - 1), 
                    burnin = 5000, iter = 5000, thin = 50)
}

# save(mixture.with.MNL.2.result, file = "mixture_with_MNL_2_result.Rdata")
# load("mixture_with_MNL_2_result.Rdata")

# convergence check
mixture.with.MNL.2.result.label <- labeling(mcmc.result = mixture.with.MNL.2.result)
for (i in 1:2) {
  print(sum(gelman.diag(mcmc.list.create(mixture.with.MNL.2.result, "beta_r", 
                                         label = mixture.with.MNL.2.result.label, 
                                         group = i), 
                        multivariate = FALSE)$psrf[, 1] > 1.1))
  print(sum(gelman.diag(mcmc.list.create(mixture.with.MNL.2.result, "beta_c", 
                                         label = mixture.with.MNL.2.result.label, 
                                         group = i), 
                        multivariate = FALSE)$psrf[, 1] > 1.1))
  print(sum(gelman.diag(mcmc.list.create(mixture.with.MNL.2.result, "sigma_r", 
                                         label = mixture.with.MNL.2.result.label, 
                                         group = i), 
                        multivariate = FALSE)$psrf[, 1] > 1.1))
  print(sum(gelman.diag(mcmc.list.create(mixture.with.MNL.2.result, "sigma_c", 
                                         label = mixture.with.MNL.2.result.label, 
                                         group = i), 
                        multivariate = FALSE)$psrf[, 1] > 1.1))
}
relabeled.gamma <- gamma.labeling(mixture.with.MNL.2.result, 
                                  mixture.with.MNL.2.result.label)
sum(gelman.diag(mcmc.list(mcmc(t(relabeled.gamma[[1]][, 1, ])), 
                          mcmc(t(relabeled.gamma[[2]][, 1, ])), 
                          mcmc(t(relabeled.gamma[[3]][, 1, ]))), 
                multivariate = FALSE)$psrf[, 1] > 1.1)

# three-group model (not identifiable)
g <- 3
mixture.with.MNL.3.result <- list()
for (i in 1:3) {
  set.seed(300 + 10 * i)
  mixture.with.MNL.3.result[[i]] <- 
    mixture_reg_MNL(y_r = y.r, y_c = y.c, x = x, z = z, 
                    id = id, l = l, g = g, 
                    b0 = matrix(0, ncol(x) + 1, g), 
                    Bvar = matrix(100, ncol(x) + 1, g), 
                    v0 = rep(0.01, g), s0 = rep(0.01, g), 
                    g0 = matrix(rep(0, (ncol(z) + 1) * (g - 1)), ncol(z) + 1, g - 1), 
                    G0 = matrix(rep(100, (ncol(z) + 1) * (g - 1)), ncol(z) + 1, g - 1), 
                    tune = matrix(rep(0.1, (ncol(z) + 1) * (g - 1)), ncol(z) + 1, g - 1), 
                    burnin = 5000, iter = 5000, thin = 50)
}

# model selection using the WAIC
WAIC(rbind(mixture.with.MNL.1.result[[1]]$loglik, 
           mixture.with.MNL.1.result[[2]]$loglik, 
           mixture.with.MNL.1.result[[3]]$loglik))
WAIC(rbind(mixture.with.MNL.2.result[[1]]$loglik, 
           mixture.with.MNL.2.result[[2]]$loglik, 
           mixture.with.MNL.2.result[[3]]$loglik))

#### postprocessing ####
# combine MCMC samples
mixture.2.beta_r.sample <- list()
for (i in 1:2) {
  mixture.2.beta_r.sample[[i]] <- t(mixture.with.MNL.2.result[[1]]$beta_r[, mixture.with.MNL.2.result.label[1, i], ])
  for (j in 2:3) {
    mixture.2.beta_r.sample[[i]] <- rbind(mixture.2.beta_r.sample[[i]], 
                                          t(mixture.with.MNL.2.result[[j]]$beta_r[, mixture.with.MNL.2.result.label[j, i], ]))
  }
}

mixture.2.beta_c.sample <- list()
for (i in 1:2) {
  mixture.2.beta_c.sample[[i]] <- t(mixture.with.MNL.2.result[[1]]$beta_c[, mixture.with.MNL.2.result.label[1, i], ])
  for (j in 2:3) {
    mixture.2.beta_c.sample[[i]] <- rbind(mixture.2.beta_c.sample[[i]], 
                                          t(mixture.with.MNL.2.result[[j]]$beta_c[, mixture.with.MNL.2.result.label[j, i], ]))
  }
}

mixture.2.sigma_r.sample <- matrix(NA, nrow(mixture.with.MNL.2.result[[1]]$sigma_r) * 3, 2)
for (i in 1:2) {
  mixture.2.sigma_r.sample[, i] <- rbind(mixture.with.MNL.2.result[[1]]$sigma_r[, mixture.with.MNL.2.result.label[1, i]], 
                                         mixture.with.MNL.2.result[[2]]$sigma_r[, mixture.with.MNL.2.result.label[2, i]], 
                                         mixture.with.MNL.2.result[[3]]$sigma_r[, mixture.with.MNL.2.result.label[3, i]])
}

mixture.2.sigma_c.sample <- matrix(NA, nrow(mixture.with.MNL.2.result[[1]]$sigma_c) * 3, 2)
for (i in 1:2) {
  mixture.2.sigma_c.sample[, i] <- rbind(mixture.with.MNL.2.result[[1]]$sigma_c[, mixture.with.MNL.2.result.label[1, i]], 
                                         mixture.with.MNL.2.result[[2]]$sigma_c[, mixture.with.MNL.2.result.label[2, i]], 
                                         mixture.with.MNL.2.result[[3]]$sigma_c[, mixture.with.MNL.2.result.label[3, i]])
}

mixture.2.G.sample <- list()
for (i in 1:2) {
  mixture.2.G.sample[[i]] <- t(mixture.with.MNL.2.result[[1]]$G[, mixture.with.MNL.2.result.label[1, i], ])
  for (j in 2:3) {
    mixture.2.G.sample[[i]] <- rbind(mixture.2.G.sample[[i]], 
                                     t(mixture.with.MNL.2.result[[j]]$G[, mixture.with.MNL.2.result.label[j, i], ]))
  }
}

# point estimates and confidence intervals of the coefficients for each group
plot.beta.summary.r <- t(rbind(apply(mixture.2.beta_r.sample[[1]], 2, 
                                     posterior.summary, CI = TRUE), 
                               apply(mixture.2.beta_r.sample[[2]], 2, 
                                     posterior.summary, CI = TRUE)))[-1, ]
plot.beta.summary.c <- t(rbind(apply(mixture.2.beta_c.sample[[1]], 2, 
                                     posterior.summary, CI = TRUE), 
                               apply(mixture.2.beta_c.sample[[2]], 2, 
                                     posterior.summary, CI = TRUE)))[-1, ]

cairo_pdf("Figure_A10.pdf", width = 6.2, height = 3.8, pointsize = 7, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0), xpd = TRUE)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location + 0.15, pch = 19, col = 1)
  points(0, plot.location - 0.15, pch = 15, col = 1)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(plot.beta.summary.c[coefficients.location, 1], plot.location + 0.15, 
           pch = 19, col = 1)
    segments(plot.beta.summary.c[coefficients.location, 2], plot.location + 0.15, 
             plot.beta.summary.c[coefficients.location, 3], plot.location + 0.15, col = 1)
    points(plot.beta.summary.c[coefficients.location, 4], plot.location - 0.15, 
           pch = 15, col = 1)
    segments(plot.beta.summary.c[coefficients.location, 5], plot.location - 0.15, 
             plot.beta.summary.c[coefficients.location, 6], plot.location - 0.15, col = 1)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend(0.4, 26, legend = c("Group 1", "Group 2"), 
       pch = c(19, 15), col = c(1, 1), bty = "n", xjust = 1, yjust = 0, cex = 0.8)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location + 0.15, pch = 19, col = 1)
  points(0, plot.location - 0.15, pch = 15, col = 1)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(plot.beta.summary.r[coefficients.location, 1], plot.location + 0.15, 
           pch = 19, col = 1)
    segments(plot.beta.summary.r[coefficients.location, 2], plot.location + 0.15, 
             plot.beta.summary.r[coefficients.location, 3], plot.location + 0.15, col = 1)
    points(plot.beta.summary.r[coefficients.location, 4], plot.location - 0.15, 
           pch = 15, col = 1)
    segments(plot.beta.summary.r[coefficients.location, 5], plot.location - 0.15, 
             plot.beta.summary.r[coefficients.location, 6], plot.location - 0.15, col = 1)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend(0.8, 26, legend = c("Group 1", "Group 2"), 
       pch = c(19, 15), col = c(1, 1), bty = "n", xjust = 1, yjust = 0, cex = 0.8)
dev.off()

# summary statistics of the intercept, error variance, 
# and proportion for each group (Table A2)
intercept.summary.c <- t(rbind(apply(mixture.2.beta_c.sample[[1]], 2, 
                                     posterior.summary, CI = TRUE), 
                               apply(mixture.2.beta_c.sample[[2]], 2, 
                                     posterior.summary, CI = TRUE)))[1, ]
round(intercept.summary.c, 2)

intercept.summary.r <- t(rbind(apply(mixture.2.beta_r.sample[[1]], 2, 
                                     posterior.summary, CI = TRUE), 
                               apply(mixture.2.beta_r.sample[[2]], 2, 
                                     posterior.summary, CI = TRUE)))[1, ]
round(intercept.summary.r, 2)

error.variance.summary.c <- apply(mixture.2.sigma_c.sample, 2, 
                                  posterior.summary, CI = TRUE)
round(error.variance.summary.c, 2)

error.variance.summary.r <- apply(mixture.2.sigma_r.sample, 2, 
                                  posterior.summary, CI = TRUE)
round(error.variance.summary.r, 2)

group.proportion.2.summary <- c(posterior.summary(rowMeans(mixture.2.G.sample[[1]]), 
                                                  CI = TRUE), 
                                posterior.summary(rowMeans(mixture.2.G.sample[[2]]), 
                                                  CI = TRUE))
round(group.proportion.2.summary, 2)

# posterior probability of individual attributions to Group 1
posterior.group.prob <- matrix(NA, ncol(mixture.2.G.sample[[1]]), 2)
for (i in 1:2) {
  posterior.group.prob[, i] <- colMeans(mixture.2.G.sample[[i]])
}

cairo_pdf("Figure_A11.pdf", width = 3.8, height = 3, pointsize = 7, family = "Helvetica")
par(mar = c(4, 4, 1, 1))
hist(posterior.group.prob[, 1], freq = FALSE, col = "white", 
     main = "", xlab = "Posterior probability to be assigned to Group 1")
dev.off()

# classification using the 5% false discovery rate threshold
group <- matrix(NA, nrow(posterior.group.prob), 2)
threshold <- rep(NA, 2)
for (i in 1:2) {
  for (j in 1:1000) {
    la <- 1 - 0.001 * j
    if (sum((1 - posterior.group.prob[, i]) * (posterior.group.prob[, i] >= la)) / 
        (sum((posterior.group.prob[, i] >= la)) + prod(posterior.group.prob[, i] < la)) > 0.05) break
  }
  la <- la + 0.001
  threshold[i] <- la
  group[, i] <- posterior.group.prob[, i] >= la
}
threshold
round(colSums(group) / nrow(posterior.group.prob), 2)
round(sum(group) / nrow(posterior.group.prob), 2)

#### postestimation simulation for the MNL ####
simulation.result <- list()
# gender
simulation.result[[1]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                          mixture.with.MNL.2.result[[2]]$gamma), 
                                        label = mixture.with.MNL.2.result.label, 
                                        data = z, colnum = 1, low = 0, high = 1, 
                                        dummy = FALSE, prob = 0.95)
# age
simulation.result[[2]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                          mixture.with.MNL.2.result[[2]]$gamma), 
                                        label = mixture.with.MNL.2.result.label, 
                                        data = z, colnum = 2, low = 2, high = 5, 
                                        dummy = FALSE, prob = 0.95)
# education
simulation.result[[3]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                          mixture.with.MNL.2.result[[2]]$gamma), 
                                        label = mixture.with.MNL.2.result.label, 
                                        data = z, colnum = 3, low = 2, high = 4, 
                                        dummy = FALSE, prob = 0.95)
# black
simulation.result[[4]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                          mixture.with.MNL.2.result[[2]]$gamma), 
                                        label = mixture.with.MNL.2.result.label, 
                                        data = z, colnum = 4, low = 0, high = 1, 
                                        dummy = TRUE, dummycol = 4:5, prob = 0.95)
# other race
simulation.result[[5]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                          mixture.with.MNL.2.result[[2]]$gamma), 
                                        label = mixture.with.MNL.2.result.label, 
                                        data = z, colnum = 5, low = 0, high = 1, 
                                        dummy = TRUE, dummycol = 4:5, prob = 0.95)
# Midwest
simulation.result[[6]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                          mixture.with.MNL.2.result[[2]]$gamma), 
                                        label = mixture.with.MNL.2.result.label, 
                                        data = z, colnum = 6, low = 0, high = 1, 
                                        dummy = TRUE, dummycol = 6:8, prob = 0.95)
# South
simulation.result[[7]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                          mixture.with.MNL.2.result[[2]]$gamma), 
                                        label = mixture.with.MNL.2.result.label, 
                                        data = z, colnum = 7, low = 0, high = 1, 
                                        dummy = TRUE, dummycol = 6:8, prob = 0.95)
# West
simulation.result[[8]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                          mixture.with.MNL.2.result[[2]]$gamma), 
                                        label = mixture.with.MNL.2.result.label, 
                                        data = z, colnum = 8, low = 0, high = 1, 
                                        dummy = TRUE, dummycol = 6:8, prob = 0.95)
# Democrat (self-identification)
simulation.result[[9]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                          mixture.with.MNL.2.result[[2]]$gamma), 
                                        label = mixture.with.MNL.2.result.label, 
                                        data = z, colnum = 9, low = 0, high = 1, 
                                        dummy = TRUE, dummycol = 9:10, prob = 0.95)
# Republican (self-identification)
simulation.result[[10]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                           mixture.with.MNL.2.result[[2]]$gamma), 
                                         label = mixture.with.MNL.2.result.label, 
                                         data = z, colnum = 10, low = 0, high = 1, 
                                         dummy = TRUE, dummycol = 9:10, prob = 0.95)
# political knowledge
simulation.result[[11]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                           mixture.with.MNL.2.result[[2]]$gamma), 
                                         label = mixture.with.MNL.2.result.label, 
                                         data = z, colnum = 11, 
                                         low = mean(z[, 11]) - sd(z[, 11]), 
                                         high = mean(z[, 11]) + sd(z[, 11]), 
                                         dummy = FALSE, prob = 0.95)
# ideology
simulation.result[[12]] <- MNL.first.dif(MCMCsample = list(mixture.with.MNL.2.result[[1]]$gamma, 
                                                           mixture.with.MNL.2.result[[2]]$gamma), 
                                         label = mixture.with.MNL.2.result.label, 
                                         data = z, colnum = 12, 
                                         low = mean(z[, 12]) - sd(z[, 12]), 
                                         high = mean(z[, 12]) + sd(z[, 12]), 
                                         dummy = FALSE, prob = 0.95)

index <- c(1:10, 12, 11)
covariate.labels <- c("Female", "Age (60s\u201330s)", 
                      "Education (BA\u2013Secondary)", 
                      "African American (Base: White)", 
                      "Other races (Base: White)", 
                      "Midwest (Base: Northeast)", 
                      "South (Base: Northeast)", 
                      "West (Base: Northeast)", 
                      "Democratic PID (Base: Indep.)", 
                      "Republican PID (Base: Indep.)", 
                      "Ideology (+1SD\u2013-1SD)", 
                      "Knowledge (+1SD\u2013-1SD)")

cairo_pdf("Figure_A12.pdf", width = 6, height = 4, pointsize = 8, family = "Helvetica")
par(mar = c(4, 0, 0, 0))
plot(NULL, NULL, bty = "n", xlim = c(-0.8, 0.4), ylim = c(0, 13),
     main = "", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 13, col = Set1[9], lwd = 0.5)
for (i in 1:12) {
  segments(-0.4, 13 - i, 0.4, 13 - i, lty = 3, col = Set1[9], lwd = 0.5)
  points(simulation.result[[index[i]]][1, 1], 13 - i, pch = 19)
  segments(simulation.result[[index[i]]][1, 2], 13 - i, 
           simulation.result[[index[i]]][1, 3], 13 - i)
  text(-0.8, 13 - i, covariate.labels[i], pos = 4)
}
mtext("First difference in probability to be assigned to Group 1", 
      side = 1, at = 0, line = 3)
axis(side = 1, at = seq(-0.4, 0.4, 0.2), lwd = 0.5)
dev.off()

# estimated differences referred in the text
round(simulation.result[[4]], 2)  # black
round(simulation.result[[10]], 2)  # Republican (self-identification)
round(simulation.result[[11]], 2)  # political knowledge
round(simulation.result[[6]], 2)  # Midwest
round(simulation.result[[8]], 2)  # West

#### results using marginal means ####
# Monte Carlo marginalization
MM.posterior.sim.result.c <- array(NA, c(20, 2, 15000))
GM.posterior.sim.result.c <- matrix(NA, 15000, 2)
for (i in 1:15000) {
  for (j in 1:2) {
    temp.data <- omit.na.data[id %in% which(mixture.2.G.sample[[j]][i, ] == 1), ]
    n.obs <- nrow(temp.data)
    GM.posterior.sim.result.c[i, j] <- mean(temp.data$choice)
    temp.mean <- tapply(temp.data$choice, temp.data$a.parents, mean)
    temp.se <- tapply(temp.data$choice, temp.data$a.parents, se)
    for (k in 1:3) {
      MM.posterior.sim.result.c[k, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
    temp.mean <- tapply(temp.data$choice, temp.data$b.identity, mean)
    temp.se <- tapply(temp.data$choice, temp.data$b.identity, se)
    for (k in 1:2) {
      MM.posterior.sim.result.c[k + 3, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
    temp.mean <- tapply(temp.data$choice, temp.data$c.activity, mean)
    temp.se <- tapply(temp.data$choice, temp.data$c.activity, se)
    for (k in 1:5) {
      MM.posterior.sim.result.c[k + 5, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
    temp.mean <- tapply(temp.data$choice, temp.data$d.policy, mean)
    temp.se <- tapply(temp.data$choice, temp.data$d.policy, se)
    for (k in 1:3) {
      MM.posterior.sim.result.c[k + 10, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
    temp.mean <- tapply(temp.data$choice, temp.data$e.habit, mean)
    temp.se <- tapply(temp.data$choice, temp.data$e.habit, se)
    for (k in 1:4) {
      MM.posterior.sim.result.c[k + 13, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
    temp.mean <- tapply(temp.data$choice, temp.data$f.intention, mean)
    temp.se <- tapply(temp.data$choice, temp.data$f.intention, se)
    for (k in 1:3) {
      MM.posterior.sim.result.c[k + 17, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
  }
  if(i %% 1000 == 0) {
    print(paste0("Iteration ", i, " is finished at ", date(), "."))
  }
}

MM.posterior.sim.result.r <- array(NA, c(20, 2, 15000))
GM.posterior.sim.result.r <- matrix(NA, 15000, 2)
for (i in 1:15000) {
  for (j in 1:2) {
    temp.data <- omit.na.data[id %in% which(mixture.2.G.sample[[j]][i, ] == 1), ]
    n.obs <- nrow(temp.data)
    GM.posterior.sim.result.r[i, j] <- mean(temp.data$rating)
    temp.mean <- tapply(temp.data$rating, temp.data$a.parents, mean)
    temp.se <- tapply(temp.data$rating, temp.data$a.parents, se)
    for (k in 1:3) {
      MM.posterior.sim.result.r[k, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
    temp.mean <- tapply(temp.data$rating, temp.data$b.identity, mean)
    temp.se <- tapply(temp.data$rating, temp.data$b.identity, se)
    for (k in 1:2) {
      MM.posterior.sim.result.r[k + 3, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
    temp.mean <- tapply(temp.data$rating, temp.data$c.activity, mean)
    temp.se <- tapply(temp.data$rating, temp.data$c.activity, se)
    for (k in 1:5) {
      MM.posterior.sim.result.r[k + 5, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
    temp.mean <- tapply(temp.data$rating, temp.data$d.policy, mean)
    temp.se <- tapply(temp.data$rating, temp.data$d.policy, se)
    for (k in 1:3) {
      MM.posterior.sim.result.r[k + 10, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
    temp.mean <- tapply(temp.data$rating, temp.data$e.habit, mean)
    temp.se <- tapply(temp.data$rating, temp.data$e.habit, se)
    for (k in 1:4) {
      MM.posterior.sim.result.r[k + 13, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
    temp.mean <- tapply(temp.data$rating, temp.data$f.intention, mean)
    temp.se <- tapply(temp.data$rating, temp.data$f.intention, se)
    for (k in 1:3) {
      MM.posterior.sim.result.r[k + 17, j, i] <- rnorm(1, temp.mean[k], temp.se[k])
    }
  }
  if(i %% 1000 == 0) {
    print(paste0("Iteration ", i, " is finished at ", date(), "."))
  }
}

# subtract the grand mean
plot.MM.summary.c <- plot.MM.summary.r <- matrix(NA, 20, 6)
for (i in 1:20) {
  plot.MM.summary.c[i, 1:3] <- posterior.summary(MM.posterior.sim.result.c[i, 1, ]) - 
    mean(GM.posterior.sim.result.c[, 1])
  plot.MM.summary.c[i, 4:6] <- posterior.summary(MM.posterior.sim.result.c[i, 2, ]) - 
    mean(GM.posterior.sim.result.c[, 2])
  plot.MM.summary.r[i, 1:3] <- posterior.summary(MM.posterior.sim.result.r[i, 1, ]) - 
    mean(GM.posterior.sim.result.r[, 1])
  plot.MM.summary.r[i, 4:6] <- posterior.summary(MM.posterior.sim.result.r[i, 2, ]) - 
    mean(GM.posterior.sim.result.r[, 2])
}

cairo_pdf("Figure_A13.pdf", width = 6.2, height = 3.8, pointsize = 7, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0), xpd = TRUE)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  for (j in 1:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(plot.MM.summary.c[coefficients.location, 1], plot.location + 0.15, 
           pch = 19, col = 1)
    segments(plot.MM.summary.c[coefficients.location, 2], plot.location + 0.15, 
             plot.MM.summary.c[coefficients.location, 3], plot.location + 0.15, col = 1, lwd = 0.8)
    points(plot.MM.summary.c[coefficients.location, 4], plot.location - 0.15, 
           pch = 15, col = 1)
    segments(plot.MM.summary.c[coefficients.location, 5], plot.location - 0.15, 
             plot.MM.summary.c[coefficients.location, 6], plot.location - 0.15, col = 1, lwd = 0.8)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend(0.4, 26, legend = c("Group 1", "Group 2"), 
       pch = c(19, 15), col = c(1, 1), bty = "n", xjust = 1, yjust = 0, cex = 0.8)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  for (j in 1:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(plot.MM.summary.r[coefficients.location, 1], plot.location + 0.15, 
           pch = 19, col = 1)
    segments(plot.MM.summary.r[coefficients.location, 2], plot.location + 0.15, 
             plot.MM.summary.r[coefficients.location, 3], plot.location + 0.15, col = 1, lwd = 0.8)
    points(plot.MM.summary.r[coefficients.location, 4], plot.location - 0.15, 
           pch = 15, col = 1)
    segments(plot.MM.summary.r[coefficients.location, 5], plot.location - 0.15, 
             plot.MM.summary.r[coefficients.location, 6], plot.location - 0.15, col = 1, lwd = 0.8)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend(0.8, 26, legend = c("Group 1", "Group 2"), 
       pch = c(19, 15), col = c(1, 1), bty = "n", xjust = 1, yjust = 0, cex = 0.8)
dev.off()

#### finite mixture linear regression models without the MNL ####
g <- 2
mixture.wo.MNL.2.result <- list()
for (i in 1:3) {
  set.seed(200 + 10 * i)
  mixture.wo.MNL.2.result[[i]] <- 
    mixture_reg(y_r = y.r, y_c = y.c, x = x, id = id, l = l, g = g, 
                b0 = matrix(0, ncol(x) + 1, g), 
                Bvar = matrix(100, ncol(x) + 1, g), 
                v0 = rep(0.01, g), s0 = rep(0.01, g), w0 = rep(1, g), 
                burnin = 5000, iter = 5000, thin = 50)
}

# save(mixture.wo.MNL.2.result, file = "mixture_wo_MNL_2_result.Rdata")
# load("mixture_wo_MNL_2_result.Rdata")

# convergence check
mixture.wo.MNL.2.result.label <- labeling(mcmc.result = mixture.wo.MNL.2.result)
for (i in 1:2) {
  print(sum(gelman.diag(mcmc.list.create(mixture.wo.MNL.2.result, "beta_r", 
                                         label = mixture.wo.MNL.2.result.label, 
                                         group = i), 
                        multivariate = FALSE)$psrf[, 1] > 1.1))
  print(sum(gelman.diag(mcmc.list.create(mixture.wo.MNL.2.result, "beta_c", 
                                         label = mixture.wo.MNL.2.result.label, 
                                         group = i), 
                        multivariate = FALSE)$psrf[, 1] > 1.1))
  print(sum(gelman.diag(mcmc.list.create(mixture.wo.MNL.2.result, "sigma_r", 
                                         label = mixture.wo.MNL.2.result.label, 
                                         group = i), 
                        multivariate = FALSE)$psrf[, 1] > 1.1))
  print(sum(gelman.diag(mcmc.list.create(mixture.wo.MNL.2.result, "sigma_c", 
                                         label = mixture.wo.MNL.2.result.label, 
                                         group = i), 
                        multivariate = FALSE)$psrf[, 1] > 1.1))
}

# postprocessing
mixture.wo.MNL.2.beta_r.sample <- list()
for (i in 1:2) {
  mixture.wo.MNL.2.beta_r.sample[[i]] <- t(mixture.wo.MNL.2.result[[1]]$beta_r[, mixture.wo.MNL.2.result.label[1, i], ])
  for (j in 2:3) {
    mixture.wo.MNL.2.beta_r.sample[[i]] <- rbind(mixture.wo.MNL.2.beta_r.sample[[i]], 
                                                 t(mixture.wo.MNL.2.result[[j]]$beta_r[, mixture.wo.MNL.2.result.label[j, i], ]))
  }
}

mixture.wo.MNL.2.beta_c.sample <- list()
for (i in 1:2) {
  mixture.wo.MNL.2.beta_c.sample[[i]] <- t(mixture.wo.MNL.2.result[[1]]$beta_c[, mixture.wo.MNL.2.result.label[1, i], ])
  for (j in 2:3) {
    mixture.wo.MNL.2.beta_c.sample[[i]] <- rbind(mixture.wo.MNL.2.beta_c.sample[[i]], 
                                                 t(mixture.wo.MNL.2.result[[j]]$beta_c[, mixture.wo.MNL.2.result.label[j, i], ]))
  }
}

mixture.wo.MNL.2.sigma_r.sample <- matrix(NA, nrow(mixture.wo.MNL.2.result[[1]]$sigma_r) * 3, 2)
for (i in 1:2) {
  mixture.wo.MNL.2.sigma_r.sample[, i] <- rbind(mixture.wo.MNL.2.result[[1]]$sigma_r[, mixture.wo.MNL.2.result.label[1, i]], 
                                                mixture.wo.MNL.2.result[[2]]$sigma_r[, mixture.wo.MNL.2.result.label[2, i]], 
                                                mixture.wo.MNL.2.result[[3]]$sigma_r[, mixture.wo.MNL.2.result.label[3, i]])
}

mixture.wo.MNL.2.sigma_c.sample <- matrix(NA, nrow(mixture.wo.MNL.2.result[[1]]$sigma_c) * 3, 2)
for (i in 1:2) {
  mixture.wo.MNL.2.sigma_c.sample[, i] <- rbind(mixture.wo.MNL.2.result[[1]]$sigma_c[, mixture.wo.MNL.2.result.label[1, i]], 
                                                mixture.wo.MNL.2.result[[2]]$sigma_c[, mixture.wo.MNL.2.result.label[2, i]], 
                                                mixture.wo.MNL.2.result[[3]]$sigma_c[, mixture.wo.MNL.2.result.label[3, i]])
}

mixture.wo.MNL.2.G.sample <- list()
for (i in 1:2) {
  mixture.wo.MNL.2.G.sample[[i]] <- t(mixture.wo.MNL.2.result[[1]]$G[, mixture.wo.MNL.2.result.label[1, i], ])
  for (j in 2:3) {
    mixture.wo.MNL.2.G.sample[[i]] <- rbind(mixture.wo.MNL.2.G.sample[[i]], 
                                            t(mixture.wo.MNL.2.result[[j]]$G[, mixture.wo.MNL.2.result.label[j, i], ]))
  }
}

# summary statistics of the intercept, error variance, 
# and proportion for each group (Table A3)
intercept.summary.c <- t(rbind(apply(mixture.wo.MNL.2.beta_c.sample[[1]], 2, 
                                     posterior.summary, CI = TRUE), 
                               apply(mixture.wo.MNL.2.beta_c.sample[[2]], 2, 
                                     posterior.summary, CI = TRUE)))[1, ]
round(intercept.summary.c, 2)

intercept.summary.r <- t(rbind(apply(mixture.wo.MNL.2.beta_r.sample[[1]], 2, 
                                     posterior.summary, CI = TRUE), 
                               apply(mixture.wo.MNL.2.beta_r.sample[[2]], 2, 
                                     posterior.summary, CI = TRUE)))[1, ]
round(intercept.summary.r, 2)

error.variance.summary.c <- apply(mixture.wo.MNL.2.sigma_c.sample, 2, 
                                  posterior.summary, CI = TRUE)
round(error.variance.summary.c, 2)

error.variance.summary.r <- apply(mixture.wo.MNL.2.sigma_r.sample, 2, 
                                  posterior.summary, CI = TRUE)
round(error.variance.summary.r, 2)

group.proportion.2.summary <- c(posterior.summary(rowMeans(mixture.wo.MNL.2.G.sample[[1]]), 
                                                  CI = TRUE), 
                                posterior.summary(rowMeans(mixture.wo.MNL.2.G.sample[[2]]), 
                                                  CI = TRUE))
round(group.proportion.2.summary, 2)

# point estimates and confidence intervals of the coefficients for each group
plot.beta.summary.r <- t(rbind(apply(mixture.wo.MNL.2.beta_r.sample[[1]], 2, 
                                     posterior.summary, CI = TRUE), 
                               apply(mixture.wo.MNL.2.beta_r.sample[[2]], 2, 
                                     posterior.summary, CI = TRUE)))[-1, ]
plot.beta.summary.c <- t(rbind(apply(mixture.wo.MNL.2.beta_c.sample[[1]], 2, 
                                     posterior.summary, CI = TRUE), 
                               apply(mixture.wo.MNL.2.beta_c.sample[[2]], 2, 
                                     posterior.summary, CI = TRUE)))[-1, ]

cairo_pdf("Figure_A14.pdf", width = 6.2, height = 3.8, pointsize = 7, family = "Helvetica")
layout(matrix(1:2, 1, 2))
par(mar = c(2, 0, 2, 0), xpd = TRUE)
plot(NULL, NULL, bty = "n", xlim = c(-1.2, 0.4), ylim = c(0, 26), 
     main = "Choice", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.4
loc.labels <- -1.07
loc.factors <- -1.2
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location + 0.15, pch = 19, col = 1)
  points(0, plot.location - 0.15, pch = 15, col = 1)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(plot.beta.summary.c[coefficients.location, 1], plot.location + 0.15, 
           pch = 19, col = 1)
    segments(plot.beta.summary.c[coefficients.location, 2], plot.location + 0.15, 
             plot.beta.summary.c[coefficients.location, 3], plot.location + 0.15, col = 1)
    points(plot.beta.summary.c[coefficients.location, 4], plot.location - 0.15, 
           pch = 15, col = 1)
    segments(plot.beta.summary.c[coefficients.location, 5], plot.location - 0.15, 
             plot.beta.summary.c[coefficients.location, 6], plot.location - 0.15, col = 1)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.4, -0.2, 0, 0.2, 0.4))
legend(0.4, 26, legend = c("Group 1", "Group 2"), 
       pch = c(19, 15), col = c(1, 1), bty = "n", xjust = 1, yjust = 0, cex = 0.8)
plot(NULL, NULL, bty = "n", xlim = c(-2.4, 0.8), ylim = c(0, 26), 
     main = "Rating", xlab = "", ylab = "", xaxt = "n", yaxt = "n")
segments(0, 0, 0, 26, col = Set1[9])
plot.location <- 26
coefficients.location <- 1
n.levels <- c(3, 2, 5, 3, 4, 3)
loc.ends <- 0.8
loc.labels <- -2.14
loc.factors <- -2.4
for (i in 1:6) {
  text(loc.factors, plot.location, labels = attribute.label[i], pos = 4)
  plot.location <- plot.location - 1
  segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
  points(0, plot.location + 0.15, pch = 19, col = 1)
  points(0, plot.location - 0.15, pch = 15, col = 1)
  text(loc.labels, plot.location, labels = level.label[[i]][1], pos = 4)
  plot.location <- plot.location - 1
  for (j in 2:n.levels[i]) {
    segments(-loc.ends, plot.location, loc.ends, plot.location, lty = 3, col = Set1[9], lwd = 0.5)
    points(plot.beta.summary.r[coefficients.location, 1], plot.location + 0.15, 
           pch = 19, col = 1)
    segments(plot.beta.summary.r[coefficients.location, 2], plot.location + 0.15, 
             plot.beta.summary.r[coefficients.location, 3], plot.location + 0.15, col = 1)
    points(plot.beta.summary.r[coefficients.location, 4], plot.location - 0.15, 
           pch = 15, col = 1)
    segments(plot.beta.summary.r[coefficients.location, 5], plot.location - 0.15, 
             plot.beta.summary.r[coefficients.location, 6], plot.location - 0.15, col = 1)
    text(loc.labels, plot.location, labels = level.label[[i]][j], pos = 4)
    plot.location <- plot.location - 1
    coefficients.location <- coefficients.location + 1
  }
}
axis(side = 1, at = c(-0.8, -0.4, 0, 0.4, 0.8))
legend(0.8, 26, legend = c("Group 1", "Group 2"), 
       pch = c(19, 15), col = c(1, 1), bty = "n", xjust = 1, yjust = 0, cex = 0.8)
dev.off()

# estimate a linear model predicting assignment to Group 2
posterior.group.prob.wo.MNL <- colMeans(mixture.wo.MNL.2.G.sample[[2]])
lm.data <- data.frame(posterior.group.prob.wo.MNL, z)
colnames(lm.data) <- c("prob.group.2", "gender", "age", "education", 
                       "black", "other.race", "Midwest", "South", "West", 
                       "Democrat", "Republican", "knowledge", "ideology")

lm.result.group.2 <- lm(prob.group.2 ~ gender + age + education + 
                          black + other.race + Midwest + South + West + 
                          Democrat + Republican + ideology + knowledge,  
                        lm.data)
round(coeftest(lm.result.group.2, vcov. = vcovHC, type = "HC3"), 3)